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Diffusion of a polymer in a gel is studied within the framework of de Gennes' model for reptation. 
Our results for the scaling of the diffusion coefficient D and the longest relaxation time r are 
markedly different from the most recently reported results, and are in agreement with de Gennes' 
reptation arguments: D ~ N~ 2 and r ~ TV 3 . The leading exponent of the finite-size corrections 
to the diffusion coefficient is consistent with the value of -2/3 that was reported for the Rubinstein 
model. This agreement suggests that its origin might be physical rather than an artifact of these 
models. 



In the theory of reptation as formulated by de GennesS 
it is predicted that the asymptotic behavior of the dif- 
fusion coefficient D of an entangled polymer scales with 
the chain length as D ~ 7V~ 2 , and that the longest re=. 
laxation time r scales as t ~ /V 3 . Evans and Edwardso 
were the first to present numerical studies of de Gennes' 
reptation model, which describes a single reptating poly- 
mer, and their results were in agreement with theoret- 
ical expectations. Recent simulations of the short-time 
behavior of this_model are also in agreement with repta- 
tion arguments]. However, the simulations of the long* 
time behavior of this model by Deutsch and Maddencl 
contradict de Gennes' predictions with, instead, D ~ 

jy-3.BiO.04 and T _ jy3.41iO.14_ 

This paper presents simulation results on the same 
model as that studied by Deutsch and Madden. The ad- 
vance in computers combined with a highly optimized 
multispin coding technique allows us to simulate over 
time scales that are several orders of magnitude longer 
than previous simulations, yielding more accurate data 
points for longer chains, which in turn allows for a data 
analysis that includes corrections to the leading power- 
law behavior. 

An alternative model for reptating polymers was pro- 
posed by RyJpjnsteinQ, and has been studied numerically 
by one of ustffl. It was found that the diffusion coefficient 
D scales as D ~ TV -2 , and the exponent of the leading 
finite-size correction was found to be close to —2/3. In an 
analytical work by Prahofer and Spohra the N~ 2 behav- 
ior was derived rigorously with a variational method, and 
the leading exponent of finite-size corrections was shown 
to be between —1/2 and —1, with the former the more 
likely. This paper will show that de Gennes' model for 
reptation and Rubinsteins repton model have the same 
length-dependence of the diffusion coefficient, both the 
leading and the subleading part in the long-length limit; 
such agreement between two different models enhances 
the credibility that these results are physical. 

In this paper, we will first describe the model consid- 
ered and then motivate the computational method used 



for the simulations. This part will be followed by the 
calculation of the diffusion constant and the results for 
the longest relaxation times. 



I. DE GENNES' MODEL FOR REPTATION 

We simulate a model for a polymer in a gel that was 
introduced by de Gennes. In this model, a polymer of 
length N is represented by N + 1 monomers on a 3D 
cubic lattice, connected by a sequence of N steps on lat- 
tice edges. One elementary move in this model is made 
by randomly selecting one monomer, and attempting to 
move it. If this monomer is an endpoint of the poly- 
mer, or if its two steps are located on the same edge, 
it will randomly move to one of the six possible lattice 
sites (the site it is already in, and five new ones). The 
time scale involved with one such an elementary move is 
At = l/(QN + 6), so that the rate of any monomer to 
move to any allowed lattice site is one time unit. 

A configuration of the polymer in the two-dimensional 
version of this model is illustrated in Figure 1. In this 
configuration, monomers 1, 5, 11, 12, and 15 can move 
to three other locations, the other monomers are frozen. 
Note that all our simulations were done in the three- 
dimensional model. 
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FIG. 1. Two-dimensional version of de Gennes' model for 
a reptating polymer with length N — 14 steps. 
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II. SIMULATION METHOD 

Rather than storing the spatial coordinates all along 
the chain, it is for simulation purposes more convenient 
to store the spatial coordinates of only the endpoints, 
and the sequence of steps from one of them to the other. 

In one elementary move in a direct implementation of 
the dynamics, the first step, the last step, or a pair of 
neighbouring steps is randomly selected and proposed to 
move. For the first and last step, this involves simply a 
replacement by a random step; a pair of steps in the mid- 
dle of the chain can only move if the steps are opposite, 
and in that case it is replaced by a randomly selected pair 
of opposite steps. Every possible change in the configu- 
ration should be proposed once per time unit, and there- 
fore 6N + 6 such elementary moves should be performed 
per unit of time. As we will see in section |y|, the typ- 
ical correlation time is approximately 0.13A^7 and thus 
around N elementary moves have to be performed to 
obtain a statistically independent configuration. Finite- 
size effects disappear only if N is well above 100, requir- 
ing arond 10 8 elementary moves per statistically indepen- 
dent configuration; an efficient implementation is crucial 
for numerical accuracy. 

To achieve such a high efficiency, we use multispin cod- 
ing. Most of our results were obtained on a DEC alpha 
workstation with a 64-bit processor, allowing 64 simu- 
lations to be run in parallel. For each of the 64 simu- 
lations, we stored the coordinates of both endpoints as 



,(°) 



,(63) 



XV' . . . X) for the left ends, and Xr . . . Xr°°' for the 
right ends. The steps ... N — 1 are stored in 3N long 
integers of 64 bits. Suppose that we denote the k th bit 
of long integer w as w k , then step i of simulation k is 
stored in the three bits {a[i] k , b[i] k , c[i] k }; if the step is in 
the positive x-, y-, or z-direction, this triplet is equal to 
{1,0,0}, {0,1,0} or {0,0,1}, respectively; the negative di- 
rections are represented by {0,1,1}, {1,0,1}, and {1,1,0}, 
respectively. Note that opposite steps are each other's 
binary complement. 

To determine whether a hairpin is located in configura- 
tion k at steps i and j = i + 1 , we have to check whether 
the triplets {a[i] k , b[i] k , c[i] k } and {a[j] fc , b[j] k , c[j] k } are 
complementary. This can be done with the following log- 
ical statement: 

h = (a[i\ 8 a[j}) A (6[i] b\j]) A (c[i] 8 c[j]); (1) 

the k th bit of long integer h is one if there is a hairpin, 
otherwise it is zero. With this one statement, we have 
however determined the location of a hairpin in all 64 
simulations! This is the power of multispin coding. 

Of course it is not sufficient to determine the presence 
of a hairpin, we also have to replace hairpins by a random 
pair of opposite steps. Suppose we have available a triplet 
of long integers {ra, rb, rc} that represent 64 random 
steps, then the logical statement 



,(63) 



and equivalent statements for b[i], c[i], a[j], b[j], and c[j] 
suffice to make one elementary move in all 64 simulations. 
In addition to these 24 logical operations some load- and 
store operations are required; altogether, we can make 64 
elementary moves (one per simulation) in a fraction of a 
microsecond. 

The endpoints require more care, since they involve 
keeping track of their spatial coordinates. For these, we 
needed a loop over the 64 simulations. 

Since we perform long simulations, the generation of 
random numbers is a major concern. In our simulation, 
there are two places where random numbers are used: 
first, in the selection of the step i, secondly in generat- 
ing a random triplet of long integers {ra, rb, rc} that 
represent 64 random steps. For the first, we used the ad- 
ditive lagged Fibonacci generator proposed by Mitchell 
and MooreO. For the second, we initially generated a set 
of M such triplets by building them up bit- by-bit with 
the Mitchell-Moore generator, and after using this set, a 
new set was generated from it by three operations: 

• each triplet was permutated with 50% probability 

• each triplet was inverted with 50 % probability 

• the list was reshuffled randomly. 

For the first and second operation, ws used a random 
bit generator proposed by Marsaglialia Note that we 
do not rely on the spectral properties of this generator, 
nor on the quality of each individual bit in the Mitchell- 
Moore generator. (Nor do we claim imperfections in these 
generators.) 



III. THE DIFFUSION CONSTANT 

To determine the diffusion constant D, we started with 
64 random polymer configurations, and let them evolve 
over a long time using the multispin approach outlined 
above. The coordinates fj and r r of the endpoints of each 
polymer were stored at regular time intervals. From their 
average f a = (fi + r r )/2 we calculate the function 

m = ((r a (tv+t)-r a (t)) 2 )-N/4 , 
6t 

the long-time limit of this function gives the diffusion 
coefficient: 



D = lim D(t). 

t — »OG 



(4) 



a[i] = (h A ra) V (^h A a[i\); 



(2) 



The additional term N/4 in Eq. (||) corrects for the rapid 
fluctuations of the endpoints relative to the center of 
mass, that do not contribute in the long-time limit; it 
reduces the effects of the finiteness of the time interval 
significantly. As we will see in section [V, the typical 
relaxation time scales as r ~ N 3 ; we find that the long- 
time limit of D{t) is well reached at t — 6A 3 , and used 



2 



this time interval to measure D. Figure 2 illustrates the 
convergence of D(t) as a function of t. 
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FIG. 2. Finite-time estimate N 2 D(t) for the rescaled dif- 
fusion constant DN 2 as a function of the length t of the time 
interval, for a polymer of length N = 60. 

The results for the diffusion coefficient D as a function 
of length TV are given in Table I. The length of the runs 
depends on the chain length N, and is also given in Table 
I. For N = 150 and N = 200, we used a 32-bit computer, 
all other results were obtained on a 64-bit computer. 

In figure ||, the measured values for DN 2 are plotted 
as a function of TV -2 / 3 ; all measurements fall on a single 
straight line given by 



DN Z = 0.173+ 1.9JV 



-2/3 



(5) 



FABLE I. The diffusion constant D times the squared 
length, for polymers of several lengths N. The third and 
fourth columns give the number of simulations and the length 
of the runs. 



N 


DN 2 


#runs 




t 


10 


0.581 


(0.001) 


G4 


3.2 


■ 10 5 


20 


0.434 


(0.001) 


64 


2.6 


■ 10 6 


30 


0.369 


(0.001) 


64 


8.6 


• 10 6 


40 


0.336 


(0.002) 


64 


4.1 


•10 7 


■50 


0.311 


(0.002) 


64 


8.0 


•10 7 


60 


0.296 


(0.003) 


64 


1.4 


• 10 s 


80 


0.275 


(0.002) 


64 


3.3 


• 10 s 


100 


0.254 


(0.005) 


64 


3.2 


• 10 s 


120 


0.251 


(0.005) 


64 


5.5 


• 10 s 


150 


0.243 


(0.004) 


32 


6.8 


• 10 s 


200 


0.226 


(0.010) 


32 


2.8 


• 10 s 
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0. 173+1. 9N~ 2/3 
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FIG. 3. Rescaled diffusion coefficient DN 2 as a function of 
jy-2/3 The g^raight line is given by equation [| 

This indicates that for long polymers, the asymptotic 
behavior of the diffusion is in agreement with reptation 
predictions of de Gennes: D ~ 7V~ 2 . Also, it gives strong 
indication that the leading exponent of the finite-size cor- 
rections has a non-integer value close to -2/3. This ex- 
ponent, ias also been observed in the Rubinstein repton 
modeB'Q, another lattice model for diffusing polymers, 
and makes it plausible that it is not simply an artifact of 
the model. 



IV. THE LONGEST RELAXATION TIMES 

The shape of the polymer at time to + 1 stays corre- 
lated to its shape at time to for a long time, but eventu- 
ally decays exponentially fast. This final decay sets the 
longest relaxation time; it is often related to the viscosity. 
The scaling of the longest relaxation time as a function 
of length is the topic of this section. To determine the 
longest relaxation time, we study the decay in the ori- 
entation d = fy — f*;. The function <f>(t) is the two-time 
autocorrelation of this vector: 



= (d(t )-d(t +t)). 



(6) 



We know that (f)(0) — N, since the expected end-to-end 
distance of the polymer is ^/N. The function <j)(i) is 
expected to decay exponentially to zero as a function of 
t, with a typical length-dependent relaxation time rjv : 



<f,(t) « ATexp(-V'rjv). 



(7) 



To determine the correlation time tjv, we have plotted in 
Figure | {-l/t')(log(<j>(t)/N - a N ) as a function of t! = 
t/N 3 ; aN is an ^-dependent parameter that removes 
partly the finite-time corrections but does not affect the 
scaling at large t': since \og((j)(t)) grows roughly linearly 
with t, it is the dominant term for longer t. The values 
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for ajv range from (X20 = —0.62 to aiso = —0.44. The 
asymptotic values of the curves in figure ^ are equal to 
N 3 /t n . 
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FIG. 4. 



For several values of N, we have plotted 



{-l/t')(log((/>(t)/N - ajv) as a function of t' = t/N 3 . All 
curves show convergence to a constant value, which is equal 
to N 3 /tjsi. From top to bottom, the curves correspond to 
N = 20, 30, 40, 50, 60, 80, 100, 120, and 150. 



It is expected from reptation arguments that in the 
limit of infinitely long chains, the product Dt scales lin- 
early with N; as we found in the previous section that 
for long chains D ~ N~ 2 , this would mean that for long 
chains t ~ N 3 . To check this, we have plotted t/N 2 as 
a function of N in figure [| For large N, the measure- 
ments are fitted well by a straight line with slope 0.13, 
indicating that for long polymers r « 0.137V 3 . 
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FIG. 5. The longest relaxation time rjv divided by iV 2 , 
as a function of TV. The straight line is a fit given by 
rjv = 0.13/V 2 +c, indicating that asymptotically r pa 0.13JV 2 . 



The values for N 3 /rjy that we obtained in this way are 
given in Table II, as well as their reciprocals; the errors 
are our estimates by eye. 



TABLE II. The longest relaxation time rjv divided by TV 3 , 
for polymers of several lengths iV. 



N 


N 3 /t n 


tn/N 3 


20 


21.2 (0.2) 


0.0472 (5) 


30 


18.0 (0.2) 


0.0556 (6) 


40 


15.7 (0.2) 


0.0637 (8) 


50 


14.2 (0.2) 


0.070 (1) 


GO 


12.9 (0.2) 


0.078 (1) 


80 


11.3 (0.2) 


0.088 (2) 


100 


10.4 (0.2) 


0.096 (2) 


120 


9.8 (0.2) 


0.102 (2) 


150 


9.3 (0.2) 


0.107 (2) 



V. CONCLUSIONS 

We have performed numerical simulations of de 
Gennes' model for polymer reptation. For the diffusion 
coefficient D, we find D = A^ 2 (0.173 + 1.9A^ 2 / 3 ). For 
the longest relaxation time r of long polymers we find 
t w 0.137V 3 . The leading-orders of the diffusion coef- 
ficient and the longest relaxation time are in agreement 
with theoretical arguments of de Gennes. There is no the- 
oretical prediction for the leading exponent of the finite- 
size corrections for the diffusion coefficient; our value of 
-2 /3 is identical to earlier reports of this exponent in the 
Rubinstein repton model. 
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